clear;clc;
global theta e12 eg ey;

Nblack = 607;
Nwhite = 3973;

rng(10010);
N_sim=100000;

wparam=[0.174027339165297;0.187646010085453;1.52703723856717;0.578390833256311;1.47397496404516;0.472254474029496;0.539039975907426;0.508222124622837;0.557092030919351;1.68356685222368;-0.456906445371444;1.53322800556343;0.520506916649071;0.554621783550490;0.495990534196635;0.550405334264551;0.503398932635956;0.149831038482197;1.21750922668894;0.774068381366255;-0.0943955231130527;-0.454412676623209;0.225964287920840;-7.47819116937451e-05;0.226994852644431;0.162646835240202];
bparam=[-0.752900651375512;-0.625236744850452;0.955164978663548;0.470520928662698;0.940898873680747;0.457615100857030;0.522441225000711;0.796501681223757;0.526383943012879;1.38074364109675;-0.825533240871413;1.01277501698395;0.503940752081312;0.230058029470066;0.270326488861762;0.435689411502657;0.141913310661526;-0.625025361512185;0.665429321866805;0.834583035216824;-0.257622993708532;-0.213896500811884;-0.532316091400443;-0.518906855782177;0.0476680272326578;0.309398571910553];

param = bparam;
b.gamma1=param(13);  b.gamma2=param(14);
b.betay=param(15);   b.cy=param(11);      b.sigtheta=param(8);    
b.c1=param(4);      b.c1d=param(21);     b.phi1=param(5);    b.phi1d=param(22);    b.beta1=param(16);  b.beta1d=param(25);
b.c2=param(9);     b.c2d=param(23);    b.phi2=param(10);   b.phi2d=param(24);   b.beta2=param(17);  b.beta2d=param(26);
b.cg=param(18);     b.cr=param(2);     b.cr=param(1);
b.phig=param(19);   b.phir=param(3);   b.phim=param(12);
b.sigg=param(20);   b.sigm=param(6);   b.sigr=param(7);


param = wparam;
w.gamma1=param(13);  w.gamma2=param(14);
w.betay=param(15);   w.cy=param(11);      w.sigtheta=param(8);    
w.c1=param(4);      w.c1d=param(21);     w.phi1=param(5);    w.phi1d=param(22);    w.beta1=param(16);  w.beta1d=param(25);
w.c2=param(9);      w.c2d=param(23);    w.phi2=param(10);   w.phi2d=param(24);   w.beta2=param(17);  w.beta2d=param(26);
w.cg=param(18);     w.cr=param(2);     w.cr=param(1);
w.phig=param(19);   w.phir=param(3);   w.phim=param(12);
w.sigg=param(20);   w.sigm=param(6);   w.sigr=param(7);


rng(10010);
ey = normrnd(0,1,N_sim,1);
theta = linspace(-2*bparam(4),2*bparam(4),50);
theta=theta';
G=b.cg+b.phig*theta;
c1tw = wparam(4)-wparam(11);
c2tw = wparam(9)-wparam(11);
c1tb = bparam(4)-bparam(11);
c2tb = bparam(9)-bparam(11);
diff=b.cy-w.cy;

% Prt1w = normcdf(bparam(4)+bparam(5)*theta+bparam(16)*G);
% Prt1w = normcdf(w.c1-w.cy+b.cy+w.phi1*theta+w.beta1*G);
Prt1w = normcdf(w.c1+diff+w.phi1*theta+w.beta1*G);
Prt1b = normcdf(w.c1+diff+w.c1d+(w.phi1+w.phi1d)*theta+(w.beta1+w.beta1d)*G);
% Prt1b = normcdf(wparam(9)+wparam(21)+(wparam(5)+wparam(22))*theta+(wparam(16)+wparam(25))*G);
figure
plot(theta,Prt1w,'k--',theta,Prt1b,'k')
xlim([-1,1]);
ylim([0,1]);
xlabel('\theta')
ylabel('Pr(T_E=1)')
legend('White Teacher, White Exp','White Teacher, Black Exp','location','southeast')
title('ELA Teacher Expectations')
saveas(gcf,'expline1.jpg')


Prt1w = normcdf(b.c1+b.phi1*theta+b.beta1*G);
Prt1b = normcdf(w.c1+diff+w.c1d+(w.phi1+w.phi1d)*theta+(w.beta1+w.beta1d)*G);
% Prt1b = normcdf(wparam(9)+wparam(21)+(wparam(5)+wparam(22))*theta+(wparam(16)+wparam(25))*G);
figure
plot(theta,Prt1w,'k--',theta,Prt1b,'k')
xlim([-1,1]);
ylim([0,1]);
xlabel('\theta')
ylabel('Pr(T_E=1)')
legend('Black Teacher, Black Exp','White Teacher, Black Exp','location','southeast')
title('ELA Teacher Expectations')
saveas(gcf,'expline3.jpg')



% Prt1w = normcdf(bparam(4)+bparam(5)*theta+bparam(16)*G);
% Prt1w = normcdf(w.c2-w.cy+b.cy+w.phi2*theta+w.beta2*G);
Prt1w = normcdf(w.c2+diff+w.phi2*theta+w.beta2*G);
Prt1b = normcdf(w.c2+diff+w.c2d+(w.phi2+w.phi2d)*theta+(w.beta2+w.beta2d)*G);
% Prt1b = normcdf(wparam(9)+wparam(21)+(wparam(5)+wparam(22))*theta+(wparam(16)+wparam(25))*G);
figure
plot(theta,Prt1w,'k--',theta,Prt1b,'k')
xlim([-1,1]);
ylim([0,1]);
xlabel('\theta')
ylabel('Pr(T_M=1)')
legend('White Teacher, White Exp','White Teacher, Black Exp','location','southeast')
title('Math Teacher Expectations')
saveas(gcf,'expline2.jpg')
% saveas(gcf,'figure1.jpg')


Prt1w = normcdf(b.c2+b.phi2*theta+b.beta2*G);
Prt1b = normcdf(w.c2+diff+w.c2d+(w.phi2+w.phi2d)*theta+(w.beta2+w.beta2d)*G);
% Prt1b = normcdf(wparam(9)+wparam(21)+(wparam(5)+wparam(22))*theta+(wparam(16)+wparam(25))*G);
figure
plot(theta,Prt1w,'k--',theta,Prt1b,'k')
xlim([-1,1])
ylim([0,1])
xlabel('\theta')
ylabel('Pr(T_M=1)')
legend('Black Teacher, Black Exp','White Teacher, Black Exp','location','southeast')
title('Math Teacher Expectations')
saveas(gcf,'expline4.jpg')









